#pragma once
#include <Dynamics/Material/PdmMaterial/pdm_linear.h>
#include <Dynamics/Material/PdmMaterial/pdm_neohookean.h>
#include <Dynamics/Material/PdmMaterial/pdm_stvk.h>
#include <Dynamics/Material/PdmMaterial/pdm_aniso.h>
#include <Dynamics/Material/PdmMaterial/pdm_nonlinear.h>
#include <Dynamics/Material/PdmMaterial/pdm_aniso_nonlinear.h>

namespace PhysLeo {

/**
 * peridynamics material type
 */
enum class PdmMaterialType
{
	LINEAR, NONLINEAR, ANISO, ANISO_NONLINEAR, STVK, NEOHOOKEAN
};

/**
 * peridynamcis material static interface.
 * \f$ W(\theta,\tau) = \frac{9\lambda}{2}g(\theta)+\frac{15\mu}{2}Gh(\tau) \f$,
 * different specific material should implement their own g,h,G function.
 */
template<typename T>
class PdmMaterial {
public:
    /**
    * return energy stored in bond according to orgin bond, deform bond, volume dilation.
    * @param[in] origin  describe the bond in rest configuration
    * @param[in] deform  describe the bond in deform configuration
    * @param[in] volume_dilation  the volume dilation of the neighborhood
    * @param[in] ptr_material  the pointer to specific vertex material data
    * @param[in] type   peridynamics material type
    * @return a T value represents the energy
    */
	__host__ __device__ static T energy(glm::tvec3<T> origin, glm::tvec3<T> deform, T volume_dilation, T* ptr_material, const PdmMaterialType type)
	{
        T a = static_cast<T>(4.5)*ptr_material[0];
        T b = static_cast<T>(7.5)*ptr_material[1];

        auto s = glm::length(deform) / glm::length(origin);
        auto si = volume_dilation;

	    switch (type)
	    {
        case PdmMaterialType::LINEAR:
        {
            typedef PdmLinear<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        case PdmMaterialType::NEOHOOKEAN:
        {
            typedef PdmNeohookean<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        case PdmMaterialType::STVK:
        {
            typedef PdmStvk<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        case PdmMaterialType::ANISO:
        {
            typedef PdmAniso<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        case PdmMaterialType::NONLINEAR:
        {
            typedef PdmNonlinear<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        case PdmMaterialType::ANISO_NONLINEAR:
        {
            typedef PdmAnisoNonlinear<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
        default:
        {
            typedef PdmLinear<T> mat;
            return a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material);
        }
	    }
	}

    /**
    * return force exerted on center point by dest point.
    * @param[in] origin  describe the bond in rest configuration
    * @param[in] deform  describe the bond in deform configuration
    * @param[in] volume_dilation  the volume dilation of the neighborhood
    * @param[in] ptr_material  the pointer to specific vertex material data
    * @param[in] type   peridynamics material type
    * @return a vec3<T> value represents the force exerted on the center point
    */
	__host__ __device__ static glm::tvec3<T> force(glm::tvec3<T> origin, glm::tvec3<T> deform, T volume_dilation, T* ptr_material, const PdmMaterialType type)
	{
        T a = static_cast<T>(4.5)*ptr_material[0];
        T b = static_cast<T>(7.5)*ptr_material[1];

        auto s = glm::length(deform) / glm::length(origin);
        auto si = volume_dilation;

        switch (type)
        {
        case PdmMaterialType::LINEAR:
        {
            typedef PdmLinear<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        case PdmMaterialType::NEOHOOKEAN:
        {
            typedef PdmNeohookean<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        case PdmMaterialType::STVK:
        {
            typedef PdmStvk<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        case PdmMaterialType::ANISO:
        {
            typedef PdmAniso<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        case PdmMaterialType::NONLINEAR:
        {
            typedef PdmNonlinear<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        case PdmMaterialType::ANISO_NONLINEAR:
        {
            typedef PdmAnisoNonlinear<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        default:
        {
            typedef PdmLinear<T> mat;
            return glm::normalize(deform)*(a * mat::dgFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::dhFunction(s, ptr_material)) / glm::length(origin);
        }
        }
	}

    /**
    * return energy stored in bond according to orgin bond, deform bond, volume dilation.
    * @param[in] origin  describe the bond in rest configuration
    * @param[in] deform  describe the bond in deform configuration
    * @param[in] volume_dilation  the volume dilation of the neighborhood
    * @param[in] ptr_material  the pointer to specific vertex material data
    * @param[in] type   peridynamics material type
    * @return a T value represents the stiffness magnitude of the bond
    */
	__host__ __device__ static T stiffnessMagnitude(glm::tvec3<T> origin, glm::tvec3<T> deform, T volume_dilation, T* ptr_material, const PdmMaterialType type)
	{
        T a = static_cast<T>(4.5)*ptr_material[0];
        T b = static_cast<T>(7.5)*ptr_material[1];

        auto s = glm::length(deform) / glm::length(origin);
        auto si = volume_dilation;

        switch (type)
        {
        case PdmMaterialType::LINEAR:
        {
            typedef PdmLinear<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        case PdmMaterialType::NEOHOOKEAN:
        {
            typedef PdmNeohookean<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        case PdmMaterialType::STVK:
        {
            typedef PdmStvk<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        case PdmMaterialType::ANISO:
        {
            typedef PdmAniso<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        case PdmMaterialType::NONLINEAR:
        {
            typedef PdmNonlinear<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        case PdmMaterialType::ANISO_NONLINEAR:
        {
            typedef PdmAnisoNonlinear<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        default:
        {
            typedef PdmLinear<T> mat;
            return (a * mat::gFunction(si, ptr_material) + b * mat::anisoWeight(glm::normalize(origin), ptr_material)*mat::hFunction(s, ptr_material)) / (glm::length(origin)*glm::length(origin));
        }
        }
	}

    /**
    * return number of material parameters.
    * @param[in] type   peridynamics material type
    * @return a int value represents number of material parameters.
    */
    __host__ __device__ static int numMaterialParams(const PdmMaterialType type)
	{
        switch (type)
        {
        case PdmMaterialType::LINEAR:
            return 2;
        case PdmMaterialType::NONLINEAR:
            return 6;
        case PdmMaterialType::ANISO:
            return 11;
        case PdmMaterialType::ANISO_NONLINEAR:
            return 15;
        case PdmMaterialType::STVK:
            return 2;
        case PdmMaterialType::NEOHOOKEAN:
            return 2;
        default:
            return 2;
        }
	}
};

}